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Abstract 

A fast numerical algorithm for the evolution of parton distributions in x 
space is described. The method is close in spirit to 'brute' force techniques. 
The necessary integrals are performed by summing the approximate contri- 
butions from small steps of the integration region. Because it is a numerical 
evaluation it shares the advantage with brute force numerical integration that 
there are no restrictions placed on the functional form of the distributions to 
be evolved. However, an improvement in the approximation technique results 
in a significant reduction in the number of integration steps and a savings in 
time on the order of three hundred fifty. The method has been implemented 
for the structure functions F% and g\ at next-to-leading order. 
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1 Introduction 



Much of the information which can be gained from studying collisions at high energy 
hadron colliders can only be extracted from the data if one has previous knowledge 
of hadronic structure. Presently this knowledge is not available by direct calculation 
and so must be obtained by other means. Dedicated high energy lepton-proton deep 
inelastic scattering, DIS, experiments provide the means to measure this structure. 
The structure functions which are measured in such experiments are essentially the 
cross-sections and hence are not directly applicable to other processes. However, they 
can be related to parton distributions which are universal and which can therefore be 
used to predict, for example, production rates in the more complicated environment 
of a hadron collider. There are other reasons why determining parton distributions 
is of interest. For example, one hopes that at some point they will be accessible by 
direct calculation. Such predictions will then need to be compared with experimental 
measurements. In addition, the proper evaluation of DIS sum rules relies on the same 
formalism by which the parton distributions are extracted. These reasons, along with 
the large and growing body of DIS data, make a strong case for efficient tools for 
determining parton distributions. 

Structure functions are related to parton distributions via the operator product 
expansion. Within this framework, the structure function at some scale, Q 2 , is ex- 
pressed in terms of a convolution of parton distributions with coefficient functions 
which can be calculated perturbatively 0. For the case of charged lepton scattering 
in the one photon exchange approximation we have, 

Fi(x, Q 2 ) = ai (x) {< e 2 > [£(Q 2 ) ® z Cf + G(Q 2 ) ® x Cf ] + q NS (Q 2 ) ® x cf s) } , 

(1) 

where / ® x q = £ ^ f {£) q{z) , a 2 = a L = x, a x = |, Z(x,Q 2 ), G(x,Q 2 ) and 
q NS (x, Q 2 ) are the quark singlet, gluon and nonsinglet distributions, respectively, the 
Ci are the coefficient functions and < e 2 > is the mean squared charge of the involved 
quarks. In ([]]) the factorization scale and the scale at which the structure function 
is defined have been implicitly set equal to each other. Because this will always be 
the case in this paper, the more general expression is not given. For g x , a± = ~ and 
the parton distributions are substituted by their polarized counterparts, AS, AG 
and Aq NS . Eq. (|I|) can be taken as a definition of the parton distributions. The 
distributions, therefore, depend on the details of the calculation of the coefficient 
functions, which are not unique. 

Though the parton distributions can not presently be calculated, their dependence 
on Q 2 is predictable for high enough values of Q 2 . Because of this, the large body of 
DIS data currently available, extending up to 5 x 10 4 GeV 2 and down to x — 10~ 5 , 
can be used together in a single global fit of the distributions. Their Q 2 dependence 
is given by the DGLAP equations @. Because the QCD processes which govern this 
dependence do not act on flavor indices, the (flavor) singlet distribution, in which these 
indices have been summed out, and the nonsinglet distribution behave differently. In 
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particular, the gluon distribution, which also does not carry flavor, can have no effect 
on the nonsinglet distribution in the perturbative evolution. Thus, in the nonsinglet 
sector we have 

^ t 1 NS M = ^P^fat)) ®. q"*®, (2) 

where t = ln^j and P NS (x, a s (Q 2 )) is the nonsinglet quark splitting function, which 
can be calculated in perturbative QCD. The coupled evolution equations in the singlet 
sector are 



dt v ; 4tt 

%,t) = ^ 
dt v ' 1 An 



P£,{<*s(t)) ®* + P qg [a s (t)) ® x G{t) 
P 9q (us{t)) ® x £(t) + P gg (a s (t)j ® x G(t) 



(3) 



The present state of knowledge of the splitting functions and coefficient functions 
allows for a next-to-leading order, NLO, treatment of both the polarized and un- 
polarized structure functions. For a concise review of the state of the art of these 
calculations see reference 0. 

The relationship (ffl-H) between the objects which are parameterized, the parton 
distributions at a fixed scale, and the data to which they are fit, structure function 
data over a wide kinematic range, is quite complicated. In fact, these equations do 
not have analytic solutions. Because of this several techniques, including integral 
transforms, polynomial expansions and numerical methods, have been employed to 
obtain approximate solutions. 

The most popular of the polynomial expansion techniques involves expanding the 
splitting and coefficient functions, as well as the parton distributions, in Laguerre 
polynomials [[§, ||. In this technique, the convolutions (|l]-||) reduce to a sum of 
products of Laguerre polynomials. A study of the Laguerre method || showed that 
one could obtain accurate results quickly in the range 0.01 < x < 0.8. Outside of this 
range, it was observed that the Laguerre expansion no longer describes the functions 
and so the method breaks down. Given that polarized scattering data is available 
down to x = 3 x 10~ 3 [7|] and unpolarized data is available down to x — 10~ 5 [^|, [|, 
this method is no longer practical. 

Another standard procedure is to transform the entire set of equations into mo- 
ment space via a Mellin transformation. The convolutions of the x dependent func- 
tions appear as products of their moment space counterparts. This simpler set of 
equations has a closed form solution which can be transformed back to x space for 
comparison with data. This technique is typically much faster than a direct numerical 
evaluation in x space, which is another popular technique. 



A 'brute force' numerical integration method was discussed in references |T(], |TT 
In this method the convolutions in Eqs. (0-0) are performed by straightforward 
numerical integration, dividing the integration region into small steps and summing 
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their approximate contributions to the integral, 

p ® x q = E r +1 -P(*)q (-) , (4) 

j Jx j Z \Z/j 

« e— ^+i^M-V (5) 

j \Xj+l/2j 

Here, P represents a general splitting or coefficient function, g is a general parton dis- 
tribution, Sxj = Xj + i — Xj and Xj+1/2 lies somewhere between and x 3+ i. Similarly, 
the differential in t is approximated by a finite difference. Eq. (0) becomes 

Ag NS (x, f J±1 ) » Ag NS (x, ft) + ^ • • E (x, +1/2 , 

(6) 

where 5tj = — £j. Similar expressions hold in the singlet sector. No restrictions 
are imposed on the functional form of the distributions because all of the functions 
are evaluated numerically. For the same reason, as was pointed out and in fact 
implemented in reference [[11], it is well suited for solving the nonlinear modified 




evolution equations jE] which take recombination effects into account. 

Because the integral on the left hand side of Eq. (f|) must be performed with 
each x step value as the lower limit, the number of times the integral on the right 
hand side must be performed is roughly proportional to the square of the number of x 
steps. This must be repeated for each point in the t grid. For the implementation of 



reference [11]] an accuracy of 2% on unpolarized distributions in the kinematic range 
relevant to the HERA collider experiments requires 1000 steps in ln(x), with a running 
time of roughly one hour.Q This is acceptable to evolve a fixed set of distributions 
but not practical for performing fits of parton distributions where the evolution may 
need to be repeated several hundred times. In addition, the freedom to try different 
functional forms, to try different data sets, or to perform other checks on the fit result 
is not easily accommodated by such a time consuming procedure. Below, a slightly 
less 'brutal' numerical technique is described which has the advantages of the brute 
force method but for which the running time is considerably reduced, making it more 
useful for performing fits. 

2 Semianalytic Solution 

Up to NLO, the splitting and coefficient functions found in the references of [0] can 
be cast in the form 

P{x) = P a (x) + £ P k fk+{x) + P S S{1 - x), (7) 

fc=0 



1 The program was run on an AlphaServer 2100 4/200. 



3 



where Po, Pi and P$ are numerical coefficients and P a (x) is a function of x. The 'plus' 
distributions, 

f ln*(l-x 

/fc+W = — i 

y 1 — sc 

are defined via the following integral over the interval to 1, 
i /ln fc (l , x , /-l ln fc (l -2 




g(z)dz= — i ifo(z) - <?(!)] dz. (9) 



\ 1 — z ) Jo I — z 

Though the method to be described does not depend on the precise form of the 
splitting and coefficient functions, it is instructive to separate the convolution integrals 
into the pieces suggested by (0). 

2.1 Basic Method and Convolutions with P„(x) 



As in the brute force method, the starting point is Eq. (Q). We would like to approx- 
imate the integrand on the right hand side of that equation as well as possible. In the 
brute force method the value of the integrand throughout the bin is approximated by 
its value at the bin center. In what will be called the semianalytic method it is first 
noted that approximating the integrals of the splitting and coefficient functions is not 
necessary because their functional dependencies are known. They could, for example, 
be explicitly integrated throughout the bin and the result multiplied by the value of 
the parton distribution at the bin center. There is a technical difficulty associated 
with this which will be described in the next section. For now we simply state that 
one solution to this difficulty is to instead approximate the distribution with a linear 
interpolation in - between its values at each pair of adjacent points in the x grid, 
z = Xj and z = Xj^\^ 



X \ H \ Xj+ , I H V x.J IX X \ IX 



+ 9 IT ■ ( 10 ) 



Substituting this expression into the right hand side of Eq. and taking only the 
P a piece of Eq. (|7J), we arrive at the expression for the convolution of P a (z) with 
q (^j in the semianalytic method, 



■ L JXn Z Jxj Z ) 



where 



Wi = ^teh^M , (12) 



Xj. 



Wo 



Xj + i Xj 
XjXj + i 



Xj + \ Xj 




(13) 
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The integrals appearing in ([11]) need only be evaluated once for each x step. The 
Wk must be evaluated for each pairwise combination of x steps at each t step of 
the evolution. If a fit is being performed, the Wk must be reevaluated each time a 
parameter is changed which affects the value of q. 

This method, in addition to retaining all of the information of the splitting and 
coefficient functions, improves the approximation of the parton distributions (which 
can only be known approximately in any case) compared to that used in the brute force 
method. The resolution of the x grid needed to obtain a reasonable approximation 
is determined by the quality of the data. Beyond a certain resolution, the accuracy 
of the approximation will exceed what the data can tell us about the distributions. 
In the implementation of this method described in which treated only g\, the 
indefinite integrals / dz^j^- and / dz^^- were explicitly worked out and were also 
evaluated numerically by gaussian quadrature. Differences between the results were 
negligible, and so in this implementation, which includes F 2 as well, all integrals 
are performed numerically. The dilogarithm, Li2(x), which occurs in the splitting 
and coefficient function has been evaluated using its power series form. Because 
performing numerical integrals involving this power series can be time consuming for 
\x\ ~ 1 the substitutions 

7T 2 

Li 2 (x) = — — ln(x)ln(l — x) — Li 2 (l — x), 

Li2(x) = _ Ll2 (0.) - Ii„* ( T l_) , 

have been used for x ~ 1 and x ~ — 1, respectively. 



2.2 Convolutions with the Plus Terms 

Apart from a minor technical detail, the plus terms are treated in the same way. 
Starting with the definition in Eq. we can write 



fk+ ®x Q 



1 , \n«(l-z) 

dz 

1-z 



1 [7, 



— q(x) / dz 



* , ln fe (l - z) 



1 - z 



(14) 



Again using Eqs. (ffl) and ( TTD ) we have 



Pkfk+® x q = Pk<J2 



W l 



ln fc (l - z) 

dz— \ + W 2 

z 1 — z) 



^■+1 , ln ft (l - z) 
dz 



z 2 1-z) 



-q[x) / dz — — 



(15) 



Because each of the integrals appearing in ([IB]) diverges if the integration region 
includes z = 1, the contribution from the last step of integration must be evaluated 
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with care. For the case of a grid of iV steps in x, the procedure is to first perform 
the sum on the right hand side of (p~5|) up to j = N — 1 and to perform the integral 
of the last term up to z = Xn only. To obtain the contribution from the N th x step, 



xjv < z < 1, the three terms on the right hand side of (|T5|) are first combined. The 
does not appear in the resulting expression, and so the integral can be explicitly 



performed. For k = and k = 1 the contributions from the N x step are 
[/o+ ®x q] N = Po [-q(x) [1 + ln{x N )) + q (^)j 



(16) 



and 



Pi [fi+ ®x q] 



N 



Pi 



x N 



I- x n 



-q(x) + q[ — 

\x N 



X 



1 — xn 



xn 



ln(l — x N ) + ln(xjv) 



+q{x) [-Li 2 (l) + U 2 {x N )} . 



(17) 



Because all of the parton distributions vanish at x = 1, their approximation with 
a linear interpolation, in addition to improving the accuracy, is responsible for the 
cancellation of the divergence. Assembling the pieces from Eqs. (|TT|), ( |i5|) and 
(|17D , and including the contribution from the delta function term in Eq. (]?]), gives 
the semianalytic expression for the convolution integrals.^ 

In , a detailed comparison of the performance of the brute force and semiana- 



lytic methods was performed for g\ in the range 0.0045 < x < 1.0 and 1 GeV < Q < 
60 GeV 2 . As the number of x steps was increased, both methods converged on the 
same result. To obtain the same accuracy within a few tenths of a percent, which is 
sufficient for the g\ data, 30 ln(Q 2 ) steps and 40 (1280) ln(x) steps were required for 
the semianalytic (brute force) methods. The semianalytic evolution was about 350 
times more efficient than the brute force method, requiring approximately 4 seconds 
to complete.^] The following section on performance will focus on the convergence of 
the semianalytic method and on the behavior of the structure functions and evolved 
parton distributions. 



3 Performance 



For the case of F 2 the MRS (A) set of distributions |T4|, defined in the MS renormaliza- 
tion scheme, was used as input to the evolution. The starting scale was 4.0 GeV 2 and 
the distributions were evaluated down to x — 10 -4 in the range 4.0 GeV 2 < Q 2 < 10 4 
GeV 2 . The evolution was performed in 200 logarithmically distributed Q 2 steps. 
Several different numbers of x steps, also distributed logarithmically, were tried. A 



2 Due to a typographical error in the presentation of this result in reference [[l3| , the full expression 
is repeated in appendix 

3 Both methods were run on a SUN SPARC 10 workstation when making this comparison. 
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reference set of distributions was first calculated using 496 x steps. Distributions 
calculated with 248, 124, and 62 x steps were then compared with these. The results 
are shown in Figs. |TJ-^for F 2 ( ), S(x) and G(x) where u v (x) is the valence up 

quark distribution and S(x) is the total sea quark distribution. 

In these figures, the difference between the reference set and the other sets has 
been normalized to the reference set and plotted as a function of x at various Q 2 
values. Note that all the quantities converge rapidly in the region from x = 10~ 4 up 
to some high value, which depends on the quantity in question. Beyond this value, 
there is a strong dependence of the result on the number of x steps. Because the 
spacing in x is logarithmic, this is not surprising. The situation could presumably be 
improved by increasing the density of points in the high x region. For example, there 
are 9 points with x > 0.5 for the case of 124 steps. Comparing the 124 step and 248 
step results for F 2 , we see that doubling this number to 18 would produce the same 
accuracy at x = 0.7 as was previously obtained at x = 0.5. This is a modest increase 
in the total number of steps and increases the running time from just under to just 
over two minutes. In || a linear step size was used at high x. 

Similar statements can be made for the parton distributions, though for the gluon 
and sea quark distributions convergence on an accurate result is not as rapid as for 
F 2 and u v . However, if one is performing a fit of parton distributions, it is only the 
accuracy with which F 2 is extracted which will influence the fit parameters. The 
convergence of the parton distributions is not directly relevant to selecting a step size 
in x and therefore does not influence the time required to perform the fit. If one only 
wants to evolve a fixed set of distributions to a high accuracy, without the iterations 
necessary to perform a fit, a large number of x steps can be used while still keeping 
a reasonable running time. In Fig. |5] a similar accuracy study is presented for g± . 

Figs. |6| and [F] show both the polarized and unpolarized parton distributions at 
two values of Q 2 . For the polarized distributions, set A from [|15| were used as 
input. As in the unpolarized case, the starting distributions were defined at 4 GeV 2 
in the MS scheme. The polarized distributions were evaluated down to x — 10~ 3 
in the range 4 GeV 2 < Q 2 < 100 GeV 2 with 62 ln(x) steps and 50 \n(Q 2 ) steps. 
The unpolarized distributions which are displayed in Fig. ^ were evaluated with 
496 ln(x) steps and 200 ln(Q 2 ) steps. The structure functions subsequently derived 
from these parton distributions are shown in Figs. |8-11|. In Figs. |^ and |9| one can 
directly see the logarithmic Q 2 behavior of the structure functions at various values 
of x. A comparison shows that the sizes of the slopes normalized to the structure 
function value are similar for the two structure functions. At the lowest values of x 
the behavior of the two structure functions is quite different; both have slopes which 
are similar in magnitude but opposite in sign. This is consistent with expectations, 
the density of partons in the low x region being dominated by the sea quark and 
gluon distributions. As the resolution of the probe increases, more partons are seen, 
but their helicity is not always remembered by the processes which generate them. 
Because the Q 2 behavior of F 2 and F\ are similar, these plots are also a dramatic 
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illustration of the danger in the low x region of assuming |r is independent of Q 2 , 
as has been done in evaluating the first moment of gi(x,Q 2 ), Ti(Q 2 ) The 
differing low x behavior of these structure functions is evident from another point 
of view in the plots of xF 2 (x) and xgi(x) in Figs. [II] and |TTJ. Here we can also see 
that the Q 2 slope of F 2 turns over around x = 0.12, which is just above the highest 
x value plotted in Fig. |[ The Q 2 slope of g\ turns over twice, behaving as F 2 in 
the valence region but becoming negative where the sea dominates. Finally, in Fig. 
[T~Q| F 2 (x,25GeV 2 ) (scaled by a factor of 50) is also plotted in order to facilitate a 



comparison with the same quantity plotted in reference 17 



4 Conclusions 

A rapid method for the evolution of parton distributions in x space has been presented. 
As a numerical technique, it places no restrictions on the form of the distributions to 
be evolved. This freedom, as well as its ability to accommodate the evolution equa- 
tions modified for recombination effects, may prove useful as the kinematic range 
of the data grows and its statistical accuracy improves. Pending a careful numer- 
ical analysis of the results, one can see now that the evolved distributions and the 
structure functions presented here compare favorably with similar figures presented 
in references Jl5| and ||17||.n 
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A Form of the Convolution Integral in the Semi- 
analytic Method 



Eq. ( |20"D gives the full expression of the semianalytic convolution of a splitting function 
or coefficient function of the form 

P(x) = P a (x) + j2Pk-( ln \ {1 _~ X) ) +P s 6(l-x), (18) 
fc=0 V 1 x / + 

with a parton distribution approximated on an x grid by 



q(-) = q<Kx /^ ■(---)+q(-), (19) 



Z J j \ Z Xa \ Xa , 

Wt-i XjJ V J/ V J/ 

4 The code may be obtained from the author via email at fasching@wisconsin.cern.ch. 
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where j is the x grid index, Xj and Xj+\ are the upper and lower boundaries of step j 
and X{ is the lower limit of the integration region. If the integration region has been 
divided into N steps in x, with x^ + i = 1, we have 



P 



x j+1 



N-1 



Psq{x 

N 

+E 

k=0 j=i 



+Pi 



Wi- 



ln fe fl - z) 



z(l-z) 



dz 



+ W 2 



ln fc (l - z) 
z 2 (l-z) 



q(xi 



, ln fc (l-^) 
1 - z 



dz 



+P \ -q( Xi ) [1 + ln(x N )\ + q I— j + q(xi)\n(l - x 

(-«{*)+ ^)) (- XN 

1 — x N V \x N /J V x 

+q(x { ) (— Li 2 (l) + Li 2 (xAr)) + ^q(xi)\n 2 (l - x<) J, 



-ln(l — x N ) + ln(xTv) 



(20) 



where Wo and are defined in Eqs. (p!2| ) and QT^). Evaluating a structure function 
at NLO at a single Q 2 value involves evaluating (p0|) for each of the values of 
for each of the involved splitting and coefficient functions. 
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Figure 1: Results for F^{x) from the evolution procedure described in the text at 
four Q 2 values. The deviation of various results from a reference result, normalized 
to the reference result, is plotted. The reference result was calculated with 496 ln(x) 
steps and 200 ln(Q 2 ) steps. The number of x steps used to obtain the other results is 
indicated on the figure. The MRS (A) |TJ] set of parton distributions at 4 GeV 2 was 
used as input. 
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Figure 2: As in Fig. 1, but for the distribution u v (x). 
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Figure 3: As in Fig. 1, but for S(x), the total sea quark distribution. 
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Figure 4: As in Fig. 1, but for G{x), the gluon distribution. 
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Figure 5: As in Fig. 1, but for the structure function gi(x). In this case the 
reference result was calculated with 248 ln(x) steps and 50 ln(Q 2 ) steps using the 



parton distributions of set A from [15] at 4 GeV as input. 
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Figure 6: Parton distributions of the proton at the starting scale of 4 GeV 2 and 
evolved to 100 GeV 2 . 
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Figure 7: The polarization of the parton distributions of the proton at the starting 
scale of 4 GeV 2 and evolved to 60 GeV 2 . 
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Figure 9: gi(Q 2 ) at four x values. 
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Figure 11: xg 1 (x) at four Q 2 values. 
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